home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / PROGRAMM / PASCAL / 0514.ZIP / CRAYZ15.ARC / VSGEDI.PAS < prev    next >
Pascal/Delphi Source File  |  1986-08-01  |  7KB  |  162 lines
  1. { Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, NY 13730 }
  2.  
  3. procedure vSGEDI ( var fID ;
  4.                     lda, n : integer ;
  5.                   var IPvt : IVECTOR ;
  6.                    var Det : DVECTOR ;
  7.                   var Work : VECTOR ;
  8.                    JobCode : integer ) ;
  9.  
  10.       { PROGRAM:                                           }
  11.       {                                                    }
  12.       {    SGEDI - computes the determinant and inverse of }
  13.       {            a matrix using the factors computed by  }
  14.       {            SGECO or SGEFA.                         }
  15.       {                                                    }
  16.       { VERSION:                 DATE:                     }
  17.       {                                                    }
  18.       {    1.4/TURBO Pascal 3.0     08/02/86               }
  19.       {    FORTRAN                  08/14/78               }
  20.       {                                                    }
  21.       { DESCRIPTION:                                       }
  22.       {                                                    }
  23.       {    On entry;                                       }
  24.       {                                                    }
  25.       {       A      - Output from SGECO or SGEFA.         }
  26.       {       lda    - Leading dimension of matrix A.      }
  27.       {       n      - Order of matrix A.                  }
  28.       {       IPvt   - Pivot index vector from SGECO or    }
  29.       {                SGEFA.                              }
  30.       {       Work   - Work vector.  Contents destroyed.   }
  31.       {       JobCode -                                    }
  32.       {            =11 compute determinant and inverse.    }
  33.       {            =01 compute inverse only.               }
  34.       {            =10 compute determinant only.           }
  35.       {                                                    }
  36.       {    On return;                                      }
  37.       {                                                    }
  38.       {       A      - Inverse of original matrix if       }
  39.       {                requested, otherwise unchanged.     }
  40.       {       Det    - Determinant of original matrix if   }
  41.       {                requested, otherwise ignored;       }
  42.       {                   determinant = Det[1]*10**Det[2]  }
  43.       {                with;                               }
  44.       {                   1.0 <= abs(Det[1]) < 10.0        }
  45.       {                or                                  }
  46.       {                   Det[1] = 0.0.                    }
  47.       {                                                    }
  48.       {    Error condition;                                }
  49.       {                                                    }
  50.       {       Division by zero will occur if the input     }
  51.       {       factor contains a zero on the diagonal and   }
  52.       {       the inverse is requested. It will not occur  }
  53.       {       if the subroutines are called correctly and  }
  54.       {       if SGECO has set RCond > 0.0 or SGEFA has    }
  55.       {       set InfoCode = 0.                            }
  56.       {                                                    }
  57.       { SUBPROGRAMS:                                       }
  58.       {                                                    }
  59.       {    SAXPY    from BLAS                              }
  60.       {    SSCAL    from BLAS                              }
  61.       {    SSWAP    from BLAS                              }
  62.       {                                                    }
  63.       { AUTHORS:                                           }
  64.       {                                                    }
  65.       {    Cleve Moler                                     }
  66.       {    University of New Mexico                        }
  67.       {    Argonne National Laboratories                   }
  68.       {                                                    }
  69.       {    PASCAL translation;                             }
  70.       {                                                    }
  71.       {    Adam Fritz                                      }
  72.       {    133 Main Street                                 }
  73.       {    Afton, New York 13730                           }
  74.  
  75. const
  76.      TEN            : real = 10.0 ;
  77.      iTEN           : integer = 10 ;
  78. var
  79.      gID            : vARRAY Absolute fID ;
  80.      d              : real ;
  81.      e              : integer ;
  82.      t              : real ;
  83.      i, j, k, l     : integer ;
  84.      kb, kp1, nm1   : integer ;
  85.  
  86. begin
  87.                                 { Compute Determinant }
  88.    if (JobCode div iTEN) <> 0 then begin
  89.       d := 1.0 ;
  90.       e := 0 ;
  91.       i := 0 ;
  92.       while (i < n) and (d <> 0.0) do begin
  93.          i := i + 1 ;
  94.          if IPvt[i] <> i then
  95.             d := -d ;
  96.          iAk := VectorRead(gID,n,i,i,1,Ak) ;
  97.          d := Ak[iAk] * d ;
  98.          if d <> 0.0 then
  99.             if Abs(d) < 1.0 then
  100.                repeat
  101.                   d := TEN * d ;
  102.                   e := e - 1
  103.                until Abs(d) >= 1.0
  104.             else if Abs(d) >= 10.0 then
  105.                repeat
  106.                   d := d/TEN ;
  107.                   e := e + 1
  108.                until Abs(d) < 10.0
  109.       end ;
  110.       Det[1] := d ;
  111.       Det[2] := e
  112.    end ;
  113.                                 { Compute inverse(U) }
  114.    if (JobCode mod iTEN) <> 0 then begin
  115.       for k := 1 to n do begin
  116.          iAk := VectorRead (gID,n,1,k,k,Ak) ;
  117.          Ak[iAk+k-1] := 1.0/Ak[iAk+k-1] ;
  118.          t := -Ak[iAk+k-1] ;
  119.          SSCAL (k-1, t, Ak[iAk], 1) ;
  120.          VectorWrite (gID,n,1,k,k,Ak) ;
  121.          if n > k then begin
  122.             kp1 := k + 1 ;
  123.             for j := kp1 to n do begin
  124.                iAj := VectorRead (gID,n,1,j,k,Aj) ;
  125.                t := Aj[iAj+k-1] ;
  126.                Aj[iAj+k-1] := 0.0 ;
  127.                SAXPY (k, t, Ak[iAk], 1, Aj[iAj], 1) ;
  128.                VectorWrite (gID,n,1,j,k,Aj) ;
  129.             end
  130.          end
  131.       end ;
  132.                                 { Form inverse(U) * inverse(L) }
  133.       if n > 1 then begin
  134.          nm1 := n - 1 ;
  135.          for kb := 1 to nm1 do begin
  136.             k := n - kb ;
  137.             kp1 := k + 1 ;
  138.             iAk := VectorRead (gID,n,1,k,n,Ak) ;
  139.             for i := kp1 to n do begin
  140.                Work[i] := Ak[iAk+i-1] ;
  141.                Ak[iAk+i-1] := 0.0
  142.             end ;
  143.             for j := kp1 to n do begin
  144.                t := Work[j] ;
  145.                iAj := VectorRead (gID,n,1,j,n,Aj) ;
  146.                SAXPY(n, t, Aj[iAj], 1, Ak[iAk], 1)
  147.             end ;
  148.             l := IPvt[k] ;
  149.             if l <> k then begin
  150.                iAl := VectorRead (gID,n,1,l,n,Al) ;
  151.                VectorWrite (gID,n,1,l,n,Ak) ;
  152.                VectorWrite (gID,n,1,k,n,Al)
  153.             end
  154.             else
  155.                VectorWrite (gID,n,1,k,n,Ak) ;
  156.          end
  157.       end
  158.    end
  159. end ;
  160.  
  161. { Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, NY 13730 }
  162.